$title  gravity.GMS


* ###################################

* Per capita income, consumption patterns, and CO2 emissions

* Journal of the Association of Environmental and Resource Economists

* Justin Caron and Thibault Fally

* April 2021

*###################################



* ESTIMATE GRAVITY EQUATIONS 

* file is capable of solving for all sectors at once for OLS
* for Poisson, solve per sector (in a loop)

$set SLASH \
$if %system.filesys% == UNIX $set SLASH /

$if not set datadir $set datadir "data%SLASH%"
* change to backslash on a PC
$setglobal datadir %datadir%

$if not set outputdatadir $set outputdatadir "estimates%SLASH%"
* change to backslash on a PC
$setglobal outputdatadir %outputdatadir%



**** MAY REQUIRE PAUSING DROPBOX


* -----------------------------------------------------------------
* DEFINE GLOBAL VARIABLES

* select dataset to use
*$if not set ds $set ds gtap7_all
*$include loaddata.gms

$if not set ds $set ds gtap8gas
$include loaddata_gtap8.gms


* GTAP9 with gas and gdt aggregated
*$if not set ds $set ds gtap9gas
*$include loaddata_gtap9.gms



* ---------------------------------------------------------------
* -- SELECT SECTORS AND REGIONS 
* ---------------------------------------------------------------

set i_(i), r_(r);
i_(i) = no; r_(r) = no;

* using all sectors :
i_(i) = yes;
*i_("nmm") = yes;


* using all countries which have distance data
r_(r)$(sum(s,importdist(r,s,"distw"))) = yes;


parameter nbr, nbi;
nbr = card(r_); nbi = card(i_);

alias(r_,s_); alias(i_,j_);

display nbr, nbi;
display r_;

parameter	coeffs parameter to store all sector specific stats
 		phiest
		logphiest
		IM
		EX
		cst;

parameter	trade(i,r,s) value of trade and absorbtion parameter in data (from to) INCLUDING INTERMEDIATE USE;

trade(i_,r_,s_) = btrade(i_,r_,s_);
trade(i_,r_,s_)$sameas(r_,s_) = sum(g,vdfm(i_,g,r_));

parameter logtrade;
logtrade(i_,r_,s_)$trade(i_,r_,s_) = log(trade(i_,r_,s_));


parameter linder;

variable sse_, loglikelihood_;


variables
delta_(*,*), IM_, EX_,
                Logfittrade(i,r,s)   defined as trade from r to s,
			tcost0(i,r,s)
			cst_(i)
			linder_(i)
			;

Equations	obj, obj_ml, tcostdef, demand,  budget, tradedef, indexdef, fmc, incomedef, cost, cost_nointer, Absorb,
		absorb_nointer, intercostdef, Cost_leontief , intercostdef_leontief, phidef ;

* NLLS : minimizing LOG error terms on trade only 
obj.. sse_ =e=      sum((i_,r_,s_)$logtrade(i_,r_,s_),  sqr( logtrade(i_,r_,s_) - logfittrade(i_,r_,s_)));

* MAXIMUM LIKELIHOOD - assuming poisson distributions

* poisson log-LIKELIHOOD function

obj_ml .. loglikelihood_ =e=  10e-3 * sum((i_,r_,s_),  trade(i_,r_,s_)*logfittrade(i_,r_,s_) - exp(logfittrade(i_,r_,s_)));


tcostdef(i_,r_,s_) .. tcost0(i_,r_,s_)  =e= 
		 -importdist(r_,s_,"ldist")*delta_(i_,"ldist")    
                +importdist(r_,s_,"contig")* delta_(i_,"contig")    
                +importdist(r_,s_,"comlang_off")* delta_(i_,"comlang_off") 
                +importdist(r_,s_,"colony")* delta_(i_,"colony")    
                +importdist(r_,s_,"homebias")* delta_(i_,"homebias")
		;

Tradedef(i_,r_,s_) ..    

	logfittrade(i_,r_,s_)
	=e= IM_(i_,s_) + EX_(i_,r_) + tcost0(i_,r_,s_) + cst_(i_) 
*	+ linder_(i_) * sqr(log(pcexp(s_)) - log(pcexp(r_))) 
;


delta_.L(i_,"ldist") = 1;		
delta_.L(i_,"contig") =	0 ;  
delta_.L(i_,"lang") = 0 ;
delta_.L(i_,"colony") =	0 ;
delta_.L(i_,"homebias") = 2;

* if using a constant :

IM_.FX(i_,"usa") = 0;
EX_.FX(i_,"usa") = 0;
Logfittrade.L(i_,r_,s_) = logtrade(i_,r_,s_);
linder_.L(i_) =0;

model calibNLP /obj, tcostdef,   tradedef/;

model calibNLP_ml /obj_ml, tcostdef,  tradedef /;


* Use PATH solver :
*option nlp = pathnlp;

* OLS
*solve calibNLP using nlp minimizing sse_;


* select set of sectors to solve :

set i__(i);


* only sectors with trade
i__(i_)$sum((r_,s_), trade(i_,r_,s_)) = yes;

*i__("tex")  = yes;
*i__("oil")  = yes;



display i__;

* -------------- -------------------------------------------------------
* solving each sector seperately :

file title_ /'title.cmd'/; title_.nd=0; title_.nw=0;

$set stime %time%


i_(i) =no;

loop(i$i__(i),

i_(i) =yes;

*       Update the title bar 
putclose title_ '@title Estimating  ',i.tl, '  (',(round(100*(ord(i)-1)/card(i))),' %% complete)  ',
		      'Start time: %stime%, Current time: %time% -- Ctrl-S to pause'/;
		      
execute 'title.cmd';


solve calibNLP_ml using nlp maximizing loglikelihood_;

coeffs(dist,"coeff",i_) = delta_.L(i_,dist) ; 

IM(i_,r_) = IM_.L(i_,r_) ;

EX(i_,r_) = EX_.L(i_,r_) ;

cst(i_) = cst_.L(i_);

linder(i_) = linder_.L(i_);

i_(i) =no;

);

parameter nbobs;

i_(i) =yes;

phiest(i_,r_) = sum(s_, exp(ex_.L(i_,s_) + tcost0.L(i_,s_,r_) + cst_.L(i_) ));
logphiest(i_,r_) = log(phiest(i_,r_));

nbobs =   sum((i_,r_,s_)$logtrade(i_,r_,s_),  1);

parameter tcostest;

tcostest(i_,s_,r_) = exp(tcost0.L(i_,s_,r_));

* check properties of poisson that sum of fitted imports and sum of true imports is equal
 parameter chkpoisson;

chkpoisson("import",r_, i_) = sum(s_, exp(logfittrade.L(i_,s_,r_)) - trade(i_,s_,r_)) ;
chkpoisson("export",r_, i_) = sum(s_, exp(logfittrade.L(i_,r_,s_)) - trade(i_,r_,s_)) ;

display chkpoisson;

parameter fittrade;

fittrade(i_,r_,s_) = exp(logfittrade.L(i_,s_,r_));

execute_unload '%outputdatadir%gravityestimates_%ds%.gdx', coeffs, im, ex,  nbobs, logphiest, phiest , cst, linder, tcostest;